# Ariga, Kenichi "Incumbency Disadvantage Under Electoral Rules with
# Intra-Party Competition: Evidence from Japan"
#
# Replication Program 1
#
# This replication program will conduct the regression discontinuity analysis at the threshold.
#
# The following figures and tables will be produced:
#      Figure 1 of the article;
#      Figures A1 and A5 through A10 of the online appendix; and
#      Tables A1 and A4 through A8 of the online appendix.

rm(list=ls())

# Set working directory.
# The outputs of this code will be saved in this directory.
setwd("")

# Load the data.
load("IncDisadvJPN.RData")

# Load the packages.
library(sandwich)
library(lmtest)

# Function to compute an Imbens-Kalyanaraman optimal bandwidth.
# Adaptation from "rdob.ado version 10.0" by Fuji, Imbens and Kalyanaraman.
# Fuji, Daisuke, Guido W. Imbens, and Karthik Kalyanaraman 2009.
# "Notes for Matlab and Stata Regression Discontinuity Software."
opt.bwd <- function(y,x,c=0){

    # Step 1: Estimation of density f(c) and conditional variance sigma-squared(c)
    # Compute the pilot bandwidth h1.
    sd.x <- sd(x)
    N <- length(x)
    h1 = 1.84*sd.x*(N^(-1/5))

    # Compute the number of units on either side of the threshold and the average outcomes.
    Nh1n <- sum( (c-h1)<=x & x<c )
    Nh1p <- sum( c<=x & x<=(c+h1) )
    Yh1n <- mean(y[(c-h1)<=x & x<c])
    Yh1p <- mean(y[c<=x & x<=(c+h1)])

    # Estimating the density of X at c and the conditional variance of Y given X_i=c at x=c
    fxc <- (Nh1p+Nh1n)/(2*N*h1)
    Y1 <- (y-Yh1n)^2
    Y1sum <- sum(Y1[(c-h1)<=x & x<c])
    Y2 <- (y-Yh1p)^2
    Y2sum <- sum(Y2[c<=x & x<=(c+h1)])
    sigma2c <- (1/(Nh1p+Nh1n))*(Y1sum+Y2sum)

    # Step 2: Estimation of second derivatives
    # Calculating the pilot bandwidth h2

    # Calculating the pilot bandwidth h2
    medXp <- median( x[x>=c] )
    medXn <- median( x[x<c] )
    Z <- data.frame(y=y,x=x,ind=(x>=c),cons=1,T=(x-c),T.sq=(x-c)^2,T.cb=(x-c)^3)
    lm1 <- lm(y ~ ind + T + T.sq + T.cb, data=Z[which(x>=medXn & x<=medXp),])
    m3c <- 6*coefficients(lm1)[5]
    Np <- sum(x>=c)
    Nn <- sum(x<c)
    h2p=3.56*((sigma2c/(fxc*max((m3c)^2, 0.01)))^(1/7))*(Np^(-1/7))
    h2n=3.56*((sigma2c/(fxc*max((m3c)^2, 0.01)))^(1/7))*(Nn^(-1/7))

    # Calculating the estimates of the second derivatives
    N2p <- sum( c<=x & x<=(c+h2p) )
    N2n <- sum( (c-h2n)<=x & x<c )
    lm.pc <- lm( y ~ T + T.sq, data=Z[which(c<=x & x<=(c+h2p)),] )
    m2pc <- 2*coefficients(lm.pc)[3]
    lm.nc <- lm( y ~ T + T.sq, data=Z[which((c-h2n)<=x & x<c),] )
    m2nc <- 2*coefficients(lm.nc)[3]

    # Step 3: Calculation of regularization terms and calculation of the optimal bandwidth

    # Calculating the regularization terms
    rp=(720*sigma2c)/(N2p*((h2p)^4))
    rn=(720*sigma2c)/(N2n*((h2n)^4))

    # Calculating the optimal bandwidth
    h.opt.unif= 5.4*(((2*sigma2c)/(fxc*(((m2pc-m2nc)^2)+(rp+rn))))^(1/5))*N^(-1/5)
    h.opt.tri= 3.4375*(((2*sigma2c)/(fxc*(((m2pc-m2nc)^2)+(rp+rn))))^(1/5))*N^(-1/5)

    h.opt <- c(h.opt.tri,h.opt.unif)
    names(h.opt) <- c("Tri","Unif")

    return(h.opt)
}

# Function to generate triangular kernel weights.
# Adaptation from "rdob.ado version 10.0" by Fuji, Imbens and Kalyanaraman.
# Fuji, Daisuke, Guido W. Imbens, and Karthik Kalyanaraman 2009.
# "Notes for Matlab and Stata Regression Discontinuity Software."
lambda <- function(x,opt.bwd) {
    z <- x/opt.bwd
    ind <- abs(z)<=1
    lambda <- (1-abs(z))*ind
    return(lambda)
}


# "party" loop begins here.
# party = 2: Conservative independents included and considered as LDP candidates.
#            "LDP" will be used in the analysis.
# party = 1: Candidates officially endorsed by LDP only. Conservative independents are excluded.
#            "LDP.oeonly" will be used in the analysis.
for (party in 1:2) {


# Load the replication data.
if (party==1) {
    master.data <- LDP.oeonly
} else if (party==2) {
    master.data <- LDP
}

# Choose the observations included in the analysis.
#
# The dataset includes elections from 1958 to 1993. The analysis requires three consecutive elections:
# t-1, t, and t+1. It estimates the causal impact of winning in t (or incumbency in t+1) on election
# outcome in t+1. Observations from t-1 are used for the placebo test.
#
# The observations for which t=1958 are deleted as they don't have t-1. Those for which t=1993 are
# also dropped as they don't have t+1.
master.data <- master.data[which(master.data$year!=1958 & master.data$year!=1993),]
# The analysis focuses on districts with 3 to 5 seats only.
master.data <- master.data[which(master.data$dm>=3 & master.data$dm<=5),]
# The observations with redistricting either between t and t-1 or between t+1 and t are dropped.
master.data <- subset(master.data, after.redist!=1)
master.data <- subset(master.data, next.after.redist!=1)
# The candidates who switched party are dropped.
master.data <- subset(master.data, next.ptyid=="LDP")
# Finally, the analysis focuses on the districts in which intra-party competition took place in t+1.
master.data <- subset(master.data, next.d.numinc>=1 & next.d.numcan.samepty > next.d.numinc)

# Set the parameters for loops.
begin.outcome <- 1   #  1 = Main analysis
end.outcome <- 3     #  2 = Outcome variables adjusted for candidates whose electoral base is inherited by an heir.
                     #  3 = Outcome variables adjusted for candidates whose electoral base is succeeded by a successor.

                     #  If outcome = 1:     If outcome = 2:          If outcome = 3:
begin.dep.var <- 1   #  1 = next.run.win    1 = next.run.win.heir    1 = next.run.win.successor
end.dep.var <- 2     #  2 = next.run.cand   2 = next.run.cand.heir   2 = next.run.cand.successor

begin.sample <- 1  # 1 = Full sample
end.sample <- 4    # 2 = Same number of candidates in elections t and t+1
                   # 3 = No suboptimal vote division in election t
                   # 4 = The above two restrictions, 2 & 3, are imposed together.

SAMPLE <- c("_a_FullSample",
            "_b_SmNC",
            "_c_NoSubopVD",
            "_d_SmNC_NoSubopVD")

# "outcome" loop begins here.
for (outcome in begin.outcome:end.outcome) {


# "sample" loop  begins here.
for (sample in 1:4) {

if (sample==1) {
    data <- master.data
} else if (sample==2) {
    data <- subset(master.data,next.d.numcan.samepty==d.numcan.samepty)
} else if (sample==3) {
    data <- subset(master.data,d.ldp.one.more.seat==0)
} else if (sample==4) {
    data <- subset(master.data,next.d.numcan.samepty==d.numcan.samepty &
                   d.ldp.one.more.seat==0)
}

# "dep.var" loop begins here.
for (dep.var in begin.dep.var:end.dep.var) {

if (outcome==1) {
    if (dep.var==1) {
        dv <- data$next.run.win
        dv.name <- "next.run.win"
    } else if (dep.var==2) {
        dv <- data$next.run.cand
        dv.name <- "next.run.cand"
    }
} else if (outcome==2) {
    if (dep.var==1) {
        dv <- data$next.run.win.heir
        dv.name <- "next.run.win.heir"
    } else if (dep.var==2) {
        dv <- data$next.run.cand.heir
        dv.name <- "next.run.cand.heir"
    }
} else if (outcome==3) {
    if (dep.var==1) {
        dv <- data$next.run.win.successor
        dv.name <- "next.run.win.successor"
    } else if (dep.var==2) {
        dv <- data$next.run.cand.successor
        dv.name <- "next.run.cand.successor"
    }
}

# Data for estimation.
est.data <- data.frame(y = dv,
                       win = data$win,
                       vt.mgn = data$d.vts.mgn,
                       candid=data$candid,
                       distid=data$distid,
                       eleid.distid=data$eleid.distid)

# I-K optimal bandwidth & triangular kernel
ik.opt.bwd <- opt.bwd(y=est.data$y,x=est.data$vt.mgn,c=0)
weights <- lambda(x=est.data$vt.mgn,ik.opt.bwd["Tri"])

# Model estimation
est.data[,"weights"] <- weights
est.data <- est.data[which(est.data$weights!=0),]
local <- lm(y ~ win + vt.mgn + win:vt.mgn,
            weights=est.data$weights, data=est.data)
coef <- coefficients(local)
N <- length(est.data$candid)
K <- local$rank
N.r <- sum(est.data$vt.mgn>0)
N.l <- sum(est.data$vt.mgn<0)

# SE clustered in candid
# Adaptation of an R code of Mahmood Arai
# http://people.su.se/~ma/mcluster.R
M <- length(unique(est.data$candid))
dfc <- (M/(M-1))*((N-1)/(N-K))
uj <- apply(estfun(local),2, function(x) tapply(x, est.data$candid, sum))
vcov.cl <- dfc*sandwich(local, meat=crossprod(uj)/N)
se.cl <- sqrt(diag(vcov.cl))
t.cl <- coef/se.cl
pvalue.cl <- 2*(1-pnorm(abs(t.cl)))
test.cl <- coeftest(local,vcov.cl)

cl <- list(vcov=vcov.cl, se=se.cl, t=t.cl, pvalue=pvalue.cl, test=test.cl, M=M)
rm(vcov.cl, se.cl, t.cl, pvalue.cl, test.cl, M, dfc, uj)

# Record the estimate of the causal effect of "win"
WIN <- c()
WIN["opt.bwd"] <- ik.opt.bwd["Tri"]
WIN["win"] <- coef["win"]
WIN["se"]     <- cl$se["win"]
WIN["t"]      <- cl$t["win"]
WIN["pvalue"] <- cl$pvalue["win"]
WIN["ub95"]   <- WIN["win"] + qnorm(0.975)*WIN["se"]
WIN["lb95"]   <- WIN["win"] + qnorm(0.025)*WIN["se"]
WIN["ub90"]   <- WIN["win"] + qnorm(0.950)*WIN["se"]
WIN["lb90"]   <- WIN["win"] + qnorm(0.050)*WIN["se"]

# Number of observations
WIN["N"] <- N
WIN["N.l"] <- N.l
WIN["N.r"] <- N.r
M1  <- length(unique(est.data$candid))
M2  <- length(unique(est.data$eleid.distid))
WIN["N.can"] <- M1
WIN["N.dist"] <- M2

if (dep.var==1) {
    RESULTS <- WIN
} else if (dep.var==2) {
    RESULTS <- rbind(RESULTS,WIN)
    if (outcome==1) {
        rownames(RESULTS) <- c("next.run.win","next.run")
    } else if (outcome==2) {
        rownames(RESULTS) <- c("next.run.win.heir","next.run.heir")
    } else if (outcome==3) {
        rownames(RESULTS) <- c("next.run.win.successor","next.run.successor")
    }
}

# "dep.var" Loop ends here.
}

RESULTS <- data.frame(RESULTS)
RESULTS$index <- c(3,2)

attach(RESULTS)

dev.new(width=10,height=2)
par(mar=c(2,7,0.25,0.25))
plot(win*100, index, type="n",
     yaxt="n", xaxt="n", xlab="", ylab="",
     ylim=c(1,4), xlim=c(-40,40)
     )
abline(v=0,lt=2,lwd=1.5,col="gray75")
points(win*100, index, pch=19, col="black")
segments(x0=lb90*100, x1=ub90*100, y0=index, y1=index,
         lty=1,lwd=2.5,col="black")
axis(1, at=seq(-40,40,10),
     labels=sprintf("%1.0f",seq(-40,40,10)),
     las=1, cex.axis=1.5)
axis(2, at=c(2,3),
     labels=c("Pr(R)","Pr(R&W)"),
     las=1, cex.axis=1.5,tick=FALSE)

if (party==2){

    if (outcome==1) {
        dev.copy(png,paste("Fig1",SAMPLE[sample],".png",sep=""),
             width=10,height=2,units="in",res=120)
        dev.off()
        if (sample==1) {
            Main.A.Full.Sample <- RESULTS
            write.csv(Main.A.Full.Sample, file = "Table_A1_a_Full_Sample.csv",row.names=TRUE)
        } else if (sample==2) {
            Main.B.SmNC <- RESULTS
            write.csv(Main.B.SmNC, file = "Table_A1_b_SmNc.csv",row.names=TRUE)
        } else if (sample==3) {
            Main.C.NoSubop.VD <- RESULTS
            write.csv(Main.C.NoSubop.VD, file = "Table_A1_c_NoSubopVD.csv",row.names=TRUE)
        } else if (sample==4) {
            Main.D.2Rest <- RESULTS
            write.csv(Main.D.2Rest, file = "Table_A1_d_2Rest.csv",row.names=TRUE)
        }
    } else if (outcome==2) {
        dev.copy(png,paste("Fig_A5_Heir",SAMPLE[sample],".png",sep=""),
             width=10,height=2,units="in",res=120)
        dev.off()
        if (sample==1) {
            Heir.A.Full.Sample <- RESULTS
            write.csv(Heir.A.Full.Sample, file = "Table_A4_Heir_a_Full_Sample.csv",row.names=TRUE)
        } else if (sample==2) {
            Heir.B.SmNC <- RESULTS
            write.csv(Heir.B.SmNC, file = "Table_A4_Heir_b_SmNc.csv",row.names=TRUE)
        } else if (sample==3) {
            Heir.C.NoSubop.VD <- RESULTS
            write.csv(Heir.C.NoSubop.VD, file = "Table_A4_Heir_c_NoSubopVD.csv",row.names=TRUE)
        } else if (sample==4) {
            Heir.D.2Rest <- RESULTS
            write.csv(Heir.D.2Rest, file = "Table_A4_Heir_d_2Rest.csv",row.names=TRUE)
        }
    } else if (outcome==3) {
        dev.copy(png,paste("Fig_A6_Suc",SAMPLE[sample],".png",sep=""),
             width=10,height=2,units="in",res=120)
        dev.off()
        if (sample==1) {
            Suc.A.Full.Sample <- RESULTS
            write.csv(Suc.A.Full.Sample, file = "Table_A5_Suc_a_Full_Sample.csv",row.names=TRUE)
        } else if (sample==2) {
            Suc.B.SmNC <- RESULTS
            write.csv(Suc.B.SmNC, file = "Table_A5_Suc_b_SmNc.csv",row.names=TRUE)
        } else if (sample==3) {
            Suc.C.NoSubop.VD <- RESULTS
            write.csv(Suc.C.NoSubop.VD, file = "Table_A5_Suc_c_NoSubopVD.csv",row.names=TRUE)
        } else if (sample==4) {
            Suc.D.2Rest <- RESULTS
            write.csv(Suc.D.2Rest, file = "Table_A5_Suc_d_2Rest.csv",row.names=TRUE)
        }
    }

} else if (party==1) {

    if (outcome==1) {
        dev.copy(png,paste("Fig_A7_OEonly",SAMPLE[sample],".png",sep=""),
             width=10,height=2,units="in",res=120)
        dev.off()
        if (sample==1) {
            OEonly.Main.A.Full.Sample <- RESULTS
            write.csv(OEonly.Main.A.Full.Sample, file = "Table_A6_OEonly_a_Full_Sample.csv",row.names=TRUE)
        } else if (sample==2) {
            OEonly.Main.B.SmNC <- RESULTS
            write.csv(OEonly.Main.B.SmNC, file = "Table_A6_OEonly_b_SmNc.csv",row.names=TRUE)
        } else if (sample==3) {
            OEonly.Main.C.NoSubop.VD <- RESULTS
            write.csv(OEonly.Main.C.NoSubop.VD, file = "Table_A6_OEonly_c_NoSubopVD.csv",row.names=TRUE)
        } else if (sample==4) {
            OEonly.Main.D.2Rest <- RESULTS
            write.csv(OEonly.Main.D.2Rest, file = "Table_A6_OEonly_d_2Rest.csv",row.names=TRUE)
        }
    } else if (outcome==2) {
        dev.copy(png,paste("Fig_A8_OEonly_Heir",SAMPLE[sample],".png",sep=""),
             width=10,height=2,units="in",res=120)
        dev.off()
        if (sample==1) {
            OEonly.Heir.A.Full.Sample <- RESULTS
            write.csv(OEonly.Heir.A.Full.Sample, file = "Table_A7_OEonly_Heir_a_Full_Sample.csv",row.names=TRUE)
        } else if (sample==2) {
            OEonly.Heir.B.SmNC <- RESULTS
            write.csv(OEonly.Heir.B.SmNC, file = "Table_A7_OEonly_Heir_b_SmNc.csv",row.names=TRUE)
        } else if (sample==3) {
            OEonly.Heir.C.NoSubop.VD <- RESULTS
            write.csv(OEonly.Heir.C.NoSubop.VD, file = "Table_A7_OEonly_Heir_c_NoSubopVD.csv",row.names=TRUE)
        } else if (sample==4) {
            OEonly.Heir.D.2Rest <- RESULTS
            write.csv(OEonly.Heir.D.2Rest, file = "Table_A7_OEonly_Heir_d_2Rest.csv",row.names=TRUE)
        }
    } else if (outcome==3) {
        dev.copy(png,paste("Fig_A9_OEonly_Suc",SAMPLE[sample],".png",sep=""),
             width=10,height=2,units="in",res=120)
        dev.off()
        if (sample==1) {
            OEonly.Suc.A.Full.Sample <- RESULTS
            write.csv(OEonly.Suc.A.Full.Sample, file = "Table_A8_OEonly_Suc_a_Full_Sample.csv",row.names=TRUE)
        } else if (sample==2) {
            OEonly.Suc.B.SmNC <- RESULTS
            write.csv(OEonly.Suc.B.SmNC, file = "Table_A8_OEonly_Suc_b_SmNc.csv",row.names=TRUE)
        } else if (sample==3) {
            OEonly.Suc.C.NoSubop.VD <- RESULTS
            write.csv(OEonly.Suc.C.NoSubop.VD, file = "Table_A8_OEonly_Suc_c_NoSubopVD.csv",row.names=TRUE)
        } else if (sample==4) {
            OEonly.Suc.D.2Rest <- RESULTS
            write.csv(OEonly.Suc.D.2Rest, file = "Table_A8_OEonly_Suc_d_2Rest.csv",row.names=TRUE)
        }
    }

}

rm(RESULTS,WIN,N,N.l,N.r,M1,M2)

# "sample" loop ends here.
}

# "outcome" loop ends here.
}


# BALANCE CHECK FOR PREDETERMINED VARIABLES
data <- master.data
begin.dep.var <- 1
end.dep.var <- 91

# "dep.var" loop begins here.
for (dep.var in begin.dep.var:end.dep.var) {

    go <- 0

    if ( (dep.var>=1 & dep.var<=9) |
        (dep.var>=21 & dep.var<=34)|
        (dep.var>=51 & dep.var<=68)|
        (dep.var>=81 & dep.var<=91)
        ) {
           go <- 1
       }

if (dep.var==1) {
    dv <- data$lag.run.cand.samedist
    dv.name <- "lag.run.cand.samedist"
} else if (dep.var==2) {
    dv <- data$lag.win.samedist
    dv.name <- "lag.win.samedist"
} else if (dep.var==3) {
    dv <- data$lag.d.vts
    dv.name <- "lag.d.vts"
} else if (dep.var==4) {
    dv <- data$numrun
    dv.name <- "numrun"
} else if (dep.var==5) {
    dv <- data$senior
    dv.name <- "senior"
} else if (dep.var==6) {
    dv <- data$age
    dv.name <- "age"
} else if (dep.var==7) {
    dv <- data$female
    dv.name <- "female"
} else if (dep.var==8) {
    dv <- data$jpmmd2010.successor
    dv.name <- "successor"
} else if (dep.var==9) {
    dv <- data$heir
    dv.name <- "heir"
}

else if (dep.var==21) {
    dv <- data$quality.cand
    dv.name <- "quality.cand"
} else if (dep.var==22) {
    dv <- data$prev.elect.office
    dv.name <- "prev.elect.office"
} else if (dep.var==23) {
    dv <- data$local.elect.office
    dv.name <- "local.elect.office"
} else if (dep.var==24) {
    dv <- data$jpmmd2010.gov
    dv.name <- "gov"
} else if (dep.var==25) {
    dv <- data$jpmmd2010.mayor
    dv.name <- "mayor"
} else if (dep.var==26) {
    dv <- data$jpmmd2010.assy
    dv.name <- "assy"
} else if (dep.var==27) {
    dv <- data$nat.loc.bur
    dv.name <- "nat.loc.bur"
} else if (dep.var==28) {
    dv <- data$jpmmd2010.bur
    dv.name <- "bur"
} else if (dep.var==29) {
    dv <- data$jpmmd2010.kenchou
    dv.name <- "kenchou"
} else if (dep.var==30) {
    dv <- data$jpmmd2010.qual
    dv.name <- "qual"
} else if (dep.var==31) {
    dv <- data$jpmmd2010.sec
    dv.name <- "sec"
} else if (dep.var==32) {
    dv <- data$jpmmd2010.news
    dv.name <- "news"
} else if (dep.var==33) {
    dv <- data$jpmmd2010.ag
    dv.name <- "ag"
} else if (dep.var==34) {
    dv <- data$jpmmd2010.talent
    dv.name <- "talent"
}

else if (dep.var==51) {
    dv <- data$dm
    dv.name <- "dm"
} else if (dep.var==52) {
    dv <- data$lag.d.ptyvts
    dv.name <- "lag.d.ptyvts"
} else if (dep.var==53) {
    dv <- data$lag.d.ptyseat.samedist
    dv.name <- "lag.d.ptyseat.samedist"
} else if (dep.var==54) {
    dv <- data$d.elect
    dv.name <- "d.elect"
} else if (dep.var==55) {
    dv <- data$lag.d.valvote
    dv.name <- "lag.d.valvote"
} else if (dep.var==56) {
    dv <- data$lag.d.turnout
    dv.name <- "lag.d.turnout"
} else if (dep.var==57) {
    dv <- data$d.numcan
    dv.name <- "d.numcan"
} else if (dep.var==58) {
    dv <- data$lag.d.eff.numcan
    dv.name <- "lag.d.eff.numcan"
} else if (dep.var==59) {
    dv <- data$d.numcan.samepty
    dv.name <- "d.numcan.samepty"
} else if (dep.var==60) {
    dv <- data$lag.d.eff.numcan.smpty.smdist
    dv.name <- "lag.d.eff.numcan.smpty.smdist"
} else if (dep.var==61) {
    dv <- data$pop.total
    dv.name <- "pop.total"
} else if (dep.var==62) {
    dv <- data$pop.female.ratio
    dv.name <- "pop.female.ratio"
} else if (dep.var==63) {
    dv <- data$pop.under15.ratio
    dv.name <- "pop.under15.ratio"
} else if (dep.var==64) {
    dv <- data$pop.over65.ratio
    dv.name <- "pop.over65.ratio"
} else if (dep.var==65) {
    dv <- data$pop.did.ratio
    dv.name <- "pop.did.ratio"
} else if (dep.var==66) {
    dv <- data$pop.work1st.ratio
    dv.name <- "pop.work1st.ratio"
} else if (dep.var==67) {
    dv <- data$pop.work2nd.ratio
    dv.name <- "pop.work2nd.ratio"
} else if (dep.var==68) {
    dv <- data$pop.work3rd.ratio
    dv.name <- "pop.work3rd.ratio"
}

else if (dep.var==81) {
    dv <- data$ele.1960
    dv.name <- "ele.1960"
} else if (dep.var==82) {
    dv <- data$ele.1963
    dv.name <- "ele.1963"
} else if (dep.var==83) {
    dv <- data$ele.1967
    dv.name <- "ele.1967"
} else if (dep.var==84) {
    dv <- data$ele.1969
    dv.name <- "ele.1969"
} else if (dep.var==85) {
    dv <- data$ele.1972
    dv.name <- "ele.1972"
} else if (dep.var==86) {
    dv <- data$ele.1976
    dv.name <- "ele.1976"
} else if (dep.var==87) {
    dv <- data$ele.1979
    dv.name <- "ele.1979"
} else if (dep.var==88) {
    dv <- data$ele.1980
    dv.name <- "ele.1980"
} else if (dep.var==89) {
    dv <- data$ele.1983
    dv.name <- "ele.1983"
} else if (dep.var==90) {
    dv <- data$ele.1986
    dv.name <- "ele.1986"
} else if (dep.var==91) {
    dv <- data$ele.1990
    dv.name <- "ele.1990"
}


# "GO" part begins here.
if (go==1) {

# Data for estimation.
est.data <- data.frame(y = dv,
                       win = data$win,
                       vt.mgn = data$d.vts.mgn,
                       candid=data$candid,
                       distid=data$distid,
                       eleid.distid=data$eleid.distid)
est.data <- na.omit(est.data)

# I-K optimal bandwidth & triangular kernel
ik.opt.bwd <- opt.bwd(y=est.data$y,x=est.data$vt.mgn,c=0)
weights <- lambda(x=est.data$vt.mgn,ik.opt.bwd["Tri"])

# Model estimation
est.data[,"weights"] <- weights
est.data <- est.data[which(est.data$weights!=0),]
local <- lm(y ~ win + vt.mgn + win:vt.mgn,
            weights=est.data$weights, data=est.data)
coef <- coefficients(local)
N <- length(est.data$candid)
K <- local$rank
N.r <- sum(est.data$vt.mgn>0)
N.l <- sum(est.data$vt.mgn<0)

M <- length(unique(est.data$candid))
dfc <- (M/(M-1))*((N-1)/(N-K))
uj <- apply(estfun(local),2, function(x) tapply(x, est.data$candid, sum))
vcov.cl <- dfc*sandwich(local, meat=crossprod(uj)/N)
se.cl <- sqrt(diag(vcov.cl))
t.cl <- coef/se.cl
pvalue.cl <- 2*(1-pnorm(abs(t.cl)))
test.cl <- coeftest(local,vcov.cl)

cl <- list(vcov=vcov.cl, se=se.cl, t=t.cl, pvalue=pvalue.cl, test=test.cl, M=M)
rm(vcov.cl, se.cl, t.cl, pvalue.cl, test.cl, M, dfc, uj)

# Record the estimate of the causal effect of "win"
WIN <- c()
WIN["opt.bwd"] <- ik.opt.bwd["Tri"]
WIN["win"] <- coef["win"]
WIN["se"]     <- cl$se["win"]
WIN["t"]      <- cl$t["win"]
WIN["pvalue"] <- cl$pvalue["win"]
WIN["ub95"]   <- WIN["win"] + qnorm(0.975)*WIN["se"]
WIN["lb95"]   <- WIN["win"] + qnorm(0.025)*WIN["se"]
WIN["ub90"]   <- WIN["win"] + qnorm(0.950)*WIN["se"]
WIN["lb90"]   <- WIN["win"] + qnorm(0.050)*WIN["se"]

# Number of observations
WIN["N"] <- N
WIN["N.l"] <- N.l
WIN["N.r"] <- N.r
M1  <- length(unique(est.data$candid))
M2  <- length(unique(est.data$eleid.distid))
WIN["N.can"] <- M1
WIN["N.dist"] <- M2

if (dep.var==1) {
    BC <- WIN
    varnames <- dv.name
} else if (dep.var>=2) {
    BC <- rbind(BC,WIN)
    varnames <- c(varnames, dv.name)
    rownames(BC) <-varnames
}

# "GO" part ends here.
}

# "dep.var" loop ends here.
}

BC <- data.frame(BC)
BC["index"] <- seq(dim(BC)[1],1,-1)

BC$pre.var <- ""

# Candidate-level variables
BC$pre.var[rownames(BC)=="lag.d.vts"] <- "Vote share in t-1"
BC$pre.var[rownames(BC)=="lag.win.samedist"] <- "Won a seat in t-1"
BC$pre.var[rownames(BC)=="senior"] <- "No. of past victories, before t"
BC$pre.var[rownames(BC)=="numrun"] <- "No. of past attempts, before t"
BC$pre.var[rownames(BC)=="age"] <- "Age"
BC$pre.var[rownames(BC)=="female"] <- "Female"
BC$pre.var[rownames(BC)=="successor"] <- 'Succeeded electoral base ("Jiban")'
BC$pre.var[rownames(BC)=="heir"] <- 'Hereditary politician'
BC$pre.var[rownames(BC)=="lag.run.cand.samedist"] <- "Ran in t-1"
BC$pre.var[rownames(BC)=="qual"] <- "Professional"
BC$pre.var[rownames(BC)=="sec"] <- "Secretary to MP"
BC$pre.var[rownames(BC)=="kenchou"] <- "Prefectural bureaucrat"
BC$pre.var[rownames(BC)=="news"] <- "Journalist"
BC$pre.var[rownames(BC)=="union"] <- "Labor union"
BC$pre.var[rownames(BC)=="ag"] <- "Agricultural coop."
BC$pre.var[rownames(BC)=="talent"] <- 'Celebrity ("Talento")'
BC$pre.var[rownames(BC)=="assy"] <- "Local representative"
BC$pre.var[rownames(BC)=="bur"] <- "National bureaucrat"
BC$pre.var[rownames(BC)=="mayor"] <- "City or municipal mayor"
BC$pre.var[rownames(BC)=="gov"] <- "Prefectural governor"
BC$pre.var[rownames(BC)=="prev.elect.office"] <- "Experience in elected office"
BC$pre.var[rownames(BC)=="local.elect.office"] <- "Locally elected official"
BC$pre.var[rownames(BC)=="quality.cand"] <- "Candidate quality"
BC$pre.var[rownames(BC)=="nat.loc.bur"] <- "National or local bureaucrat"

# District-level variables
BC$pre.var[rownames(BC)=="dm"] <- "District Magnitude"
BC$pre.var[rownames(BC)=="lag.d.ptyvts"] <- "LDP's vote share in t-1"
BC$pre.var[rownames(BC)=="lag.d.ptyseat.samedist"] <- "No. of seats won by LDP in t-1"
BC$pre.var[rownames(BC)=="d.elect"] <- "Size of electorate"
BC$pre.var[rownames(BC)=="d.valvote"] <- "Total votes cast in t"
BC$pre.var[rownames(BC)=="d.turnout"] <- "Turnout in t"
BC$pre.var[rownames(BC)=="lag.d.valvote"] <- "Total votes cast in t-1"
BC$pre.var[rownames(BC)=="lag.d.turnout"] <- "Turnout in t-1"
BC$pre.var[rownames(BC)=="d.numcan"] <- "No. of candidates in t"
BC$pre.var[rownames(BC)=="d.eff.numcan"] <- "Eff. no. of candidates in t"
BC$pre.var[rownames(BC)=="lag.d.eff.numcan"] <- "Eff. no. of candidates in t-1"
BC$pre.var[rownames(BC)=="d.numcan.samepty"] <- "No. of LDP candidates in t"
BC$pre.var[rownames(BC)=="d.eff.numcan.samepty"] <- "Eff. no. of LDP candidates in t"
BC$pre.var[rownames(BC)=="lag.d.eff.numcan.smpty.smdist"] <- "Eff. no. of LDP candidates in t-1"
BC$pre.var[rownames(BC)=="d.numpty"] <- "No. of parties in t"
BC$pre.var[rownames(BC)=="d.eff.numpty"] <- "Eff. no. of parties in t"
BC$pre.var[rownames(BC)=="lag.d.eff.numpty"] <- "Eff. no. of parties in t-1"
BC$pre.var[rownames(BC)=="pop.total"] <- "District population"
BC$pre.var[rownames(BC)=="pop.male.ratio"] <- "Proportion of male residents"
BC$pre.var[rownames(BC)=="pop.female.ratio"] <- "Proportion of female residents"
BC$pre.var[rownames(BC)=="pop.under15.ratio"] <- "Residents under 15 years old"
BC$pre.var[rownames(BC)=="pop.over65.ratio"] <- "Residents over 65 years old"
BC$pre.var[rownames(BC)=="pop.did.ratio"] <- "Proportion of urban residents"
BC$pre.var[rownames(BC)=="pop.work1st.ratio"] <- "Primary sector employment"
BC$pre.var[rownames(BC)=="pop.work2nd.ratio"] <- "Industry sector employment"
BC$pre.var[rownames(BC)=="pop.work3rd.ratio"] <- "Service sector employment"

# Election year dummies
BC$pre.var[rownames(BC)=="ele.1960"] <- "Election t = 1960"
BC$pre.var[rownames(BC)=="ele.1963"] <- "Election t = 1963"
BC$pre.var[rownames(BC)=="ele.1967"] <- "Election t = 1967"
BC$pre.var[rownames(BC)=="ele.1969"] <- "Election t = 1969"
BC$pre.var[rownames(BC)=="ele.1972"] <- "Election t = 1972"
BC$pre.var[rownames(BC)=="ele.1976"] <- "Election t = 1976"
BC$pre.var[rownames(BC)=="ele.1979"] <- "Election t = 1979"
BC$pre.var[rownames(BC)=="ele.1980"] <- "Election t = 1980"
BC$pre.var[rownames(BC)=="ele.1983"] <- "Election t = 1983"
BC$pre.var[rownames(BC)=="ele.1986"] <- "Election t = 1986"
BC$pre.var[rownames(BC)=="ele.1990"] <- "Election t = 1990"

dev.new(width=5,height=10)
par(mar=c(5,9,1,1))
    plot(BC$pvalue,BC$index, type="n",
         yaxt="n", xaxt="n", xlab="p-value", ylab="",
         xlim=c(0,1))
for (i in 1:length(BC$pvalue)) {
    abline(h=i, lty=1, col="gray75")
}
    points(BC$pvalue,BC$index, pch=19, cex=0.5, col="black")
abline(v=0.05, lty=2, col="black")
abline(v=0.10, lty=3, col="black")
axis(1, at=c(seq(0,0.1,0.05),seq(0.1,1,0.1)),
     labels=sprintf("%1.2f",c(seq(0,0.1,0.05),seq(0.1,1,0.1))), las=2, cex.axis=0.75)
axis(2, at=BC$index, labels=BC$pre.var, las=2, tick=FALSE, cex.axis=0.6)

if (party==2) {
    dev.copy(png,"Fig_A1_PlaceboTest.png",
                 width=5,height=10,units="in",res=120)
} else if (party==1){
    dev.copy(png,"Fig_A10_PlaceboTest_OEonly.png",
                 width=5,height=10,units="in",res=120)
}
dev.off()


graphics.off()

# "party" loop ends here.
}
